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1.  INTRODUCTION 

Figure  1  shows  two  output0  from  the  ray  tracing  subroutine  of  the  common  shot 
inversion  code,  CXZCS  (Dong,  1989).  In  both  diagrams,  the  ray  from  the  source  to 
the  output  point  and  a  suite  of  rays  from  the  output  point  to  a  suite  of  receivers  are 
shown.  The  spread  of  rays  in  the  lower  figure  (Figure  lb)  can  be  seen  to  be  smoother 
than  the  spread  in  the  upper  figure  (Figure  la).  While  the  travel  times  for  the  two 
sets  of  rays  is  nearly  the  same,  the  amplitude-versus-offset  distribution  for  the  rays 
of  Figure  lb  is  also  smoother  than  the  amplitude  distribution  for  the  rays  of  Figure 
la. 

The  model  for  the  Figure  la  was  generated  by  picking  points  on  the  interfaces 
and  then  connecting  them  with  cubic  splines.  For  Figure  lb,  the  model  of  Figure 
la  is  passed  through  a  smoothing  operator  which  penalizes  point  distributions  that 
produce  large  derivatives,  while  still  requires  that  the  modified  curve  describing  each 
interface  still  remains  "close”  to  the  original  model.  This  paper  describes  the  smooth¬ 
ing  operator  used  to  generate  the  model  used  in  Figure  lb. 

The  model  generating  techinique  described  above  is  an  example  of  the  process  of 
observing  and  recording  data — in  that  case,  (x,  z)  points  of  model  interfaces.  Usually, 
these  data  are  inaccurate  because  of  noise  or  measurement  error.  If  we  directly  connect 
these  data  by  interpolation  to  obtain  a  curve,  the  result  often  exhibits  a  roughness 
that  characterizes  the  inaccuracies  due  to  the  errors,  and  suggests  a  fineness  of  detail 
that  is  invalid  for  the  data.  When  one  has  a  priori  information  that  data  represents 
a  smooth  phenomenon,  then  it  is  important  to  smooth  the  observed  data  consistent 
with  the  priori  knowlwdge. 

Two  common  smoothing  methods  are:  (1)  lowpass  Fourier  transform  and  (2)  the 
Moving-average  Method.  Lowpass  Fourier  transform  is  carried  out  in  the  frequency 
domain.  This  method  removes  all  frequencies  above  a  specific  frequency.  We  can  use 
FFT  to  speed  up  the  computation,  but  it  requires  uniform  grid  data.  The  Moving- 
average  Method  replaces  a  value  at  a  point  by  an  average  value  around  that  point,  but 
it  cannot  suppress  high  frequency  components  efficiently.  In  this  paper,  we  give  an 
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Fig.  1.  Ray  tracing  in  the  common  shot  data.  The  lower  one  is  the  smoothed  model. 
The  non-smooth  model  has  one  ray  tracing  broken. 


2 


Liu 


Smoothing  Technique 


alternative  smoothing  method.  A  smooth  function  must  have  small  derivatives.  So  we 
hope  that  a  smooth  curve  is  both  close  to  the  original  one  and  has  as  small  derivatives 
as  possible.  That  is,  we  minimize  the  weighted  squared  sum  of  the  error  between  the 
desired  curve  and  the  original  one  and  some  order  derivative  of  the  desired  curve. 
This  method  suppresses  high  frequency  components  efficiently  instead  of  removing 
them.  In  discretization,  this  method  does  not  require  uniform  grid  data,  and  the 
smooth  curve  is  found  by  solving  a  banded  linear  system  whose  computational  cost 
is  proportional  to  number  of  data  points. 

2.  SMOOTHING  TECHNIQUE 
Mathematical  Equation 

Let  z  —  f(x)  be  any  continuous  function.  We  solve  for  a  smooth  function  g(x) 
that  approximates  f(x)  through  the  requirement 

f[f(x)-g(*)]2dx  +  a min,  (1) 

where  a  >  0  is  called  the  smoothing  parameter  and  n  is  a  positive  integer.  The  larger 
the  value  of  a  and  the  larger  n,  the  smoother  g(x)  will  be. 

By  using  the  calculus  of  variations  we  can  change  (1)  into  a  differential  equation 
for  g(x): 

(-')"»  0  +  pW  =  /W. 

This  equation  has  a  simple  solution  by  Fourier  transform: 

G(k)  =  F(*)/(  1  +  (2) 

Here,  we  have  used  capitol  letters  to  denote  the  Fourier  transforms  of  the  correspond¬ 
ing  functions  denoted  by  lower  case  letters  in  the  spatial  domain.  This  expression 
shows  that  high- wavenumber  components  of  f(x)  are  suppressed  in  the  approximation 
g{x).  Equation  (2)  is  same  as  the  Butterworth  filter. 

Numerical  Algorithm 

In  general,  we  know  a  discrete  data  set,  (x,,  /,),  i  =  0, 1,2,...,  N,  instead  of  con¬ 
tinuous  data,  f(x).  Let  /t*  =  x<  —  x*_i.  We  use  a  finite  difference  approximation 
of  the  differential  operator  in  equation  (1).  If  the  h,-’s  are  equal,  this  work  is  not 
difficult.  Otherwise,  it  becomes  very  complicated  for  n  more  than  2.  Fortunately, 
n  =  2  is  a  smooth  enough  index.  After  discretization  equation  (!)  is  equivalent  to 
a  linear  algebric  system  of  equations  for  the  unknown  g,'s.  The  coefficient  matrix  of 
the  system  is  symmetric  positive  definite,  for  each  n.  When  n  =  1  the  matrix  is 
tridiagonal;  when  n  =  2,  the  matrix  is  five-diagonal.  Their  comp'ffations  for  solving 
the  g,'s  are  proportional  to  N. 
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Choice  of  smoothing  parameter 

If  the  smoothing  parameter  a  is  two  large,  we  will  lose  some  low  frequency  compo¬ 
nents.  Therefore,  choosing  a  suitable  smoothing  parameter  is  necessary.  We  specify 
a  wavenumber  k0,  then  choose  a  such  that 

G(*o)  =  0.5  F(k0). 

Here  k0  should  be  chosen  as  some  proportion  of  the  Nyquist  wavenumber,  kNyq.  For 
example,  empirically,  I  found  that  taking  ko  =  A'jvV9/6,  produced  travel  times  not 
too  far  from  the  traveltimes  of  the  original  model,  but  amplitudes  that  were  much 
smoother  as  a  function  of  receiver  position  than  for  the  unsmoothed  model  interfaces. 


3.  COMPARISON  OF  SMOOTHING  METHODS 


The  Moving-average  Method  smooths  a  function  as 


9(x) 


1 

2&x 


(3) 


where  f(x)  is  an  original  function  and  p(x)  is  the  desired  function.  In  the  wavenumber 
domain,  the  solution  is 


G(k)  = 


sin(A:A:r) 

kAx 


F(k). 


The  larger  Ax  is,  the  smoother  g(x)  will  be. 


(4) 


Fig.  2.  Smoothing  filters  in  the  wavenumber  domain.  Label  M  represents  the 
Moving- average  method;  n  =  1,2  represent  our  smoothing  method  with  n  as  defined 
in  equation  (1).  L  represents  the  lowpass  Fourier  transform. 

From  Figure  2,  lowpass  Fourier  transform  can  be  seen  to  most  effectively  eliminate 
high  frequencies,  however,  at  a  cost  of  producing  ringing  in  the  spatial  domain.  With 


4 


Liu 


Smoothing  Technique 


our  method,  n  =  2  provides  a  sufficiently  good  approximation  of  the  rectangle  of 
the  lowpass  Fourier  transform  while  not  producing  the  same  undesirable  effect  in  the 
spatial  domain.  The  (three-term)  Moving-average  filter  decays  most  slowly  in  the 
wavenumber  domain.  This  filter  with  larger  windows  and  nonuniform  sampling  is 
more  cunbersome  to  use  than  the  filter  we  propose  here. 

In  computational  cost,  Fourier  transform  is  0(N  log  N),  the  Moving-average  is 
O(N) — but  with  a  large  constant  when  decay  rate  in  the  Fourier  domain  is  made  to 
compete  with  our  approach,  and  our  method  is  0(n2N). 

4.  APPLICATION  TO  RAY  TRACING 

In  seismic  data  processing,  we  often  meet  ray  tracing  problems.  In  particular,  we 
consider  the  medium  to  consist  of  constant-velocity  layers  separated  by  arbitrary  in¬ 
terfaces  which  are  obtained  by  measurement  or  digitizing  from  graphs.  To  guarantee 
a  successful  ray  tracing,  the  modeled  interface  must  be  made  smooth.  Two  examples 
are  illustrated  in  Figure  1  and  3.  The  smoothed  interfaces  yield  a  more  uniform  cover¬ 
age  by  rays,  which  gurantees  a  stable  amplitude  computation  (Liu,  1991).  Moreover, 
the  smoothed  interfaces  avoid  the  ray  tracing  broken  (Figure  1). 

5.  CONCLUSIONS 

In  this  paper,  I  introduced  a  smoothing  technique  for  experimental  curve.  This 
method  is  not  restricted  by  uniform  sampling,  and  smooths  a  curve  efficiently.  A  key 
step  in  practice  is  how  to  choose  a  suitable  smoothing  parameter.  This  depends  on 
the  properties  of  the  problem  to  be  solved.  In  the  example  of  my  application,  I  chose 
the  smoothing  parameter  such  that  the  travaltime  is  preserved  and  the  amplitude 
ditribution  becomes  smoother. 
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